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ABSTRACT 

We present a merger-driven evolutionary model for the production of luminous (LIRGs) 
and ultraluminous infrared galaxies (ULIRGs) with warm IR colours. Our results show that 
simulations of gas-rich major mergers including star formation, black hole growth, and feed- 
back can produce warm (U)LIRGs. We also find that while the warm evolutionary phase is 
associated with increased AGN activity, star formation alone may be sufficient to produce 
warm IR colours. However, the transition can be suppressed entirely - even when there is a 
significant AGN contribution - when we assume a single-phase ISM, which maximizes the 
attenuation. Finally, our evolutionary models are consistent with the 25-to-60 ftm flux density 
ratio vs. Lhx /Ljr relation for local LIRGs and ULIRGs, and predict the observed scatter 
in IR color at fixed Lhx /Liu. Therefore, our models suggest a cautionary note in the inter- 
pretation of warm IR colours: while associated with periods of active black hole growth, they 
are probably produced by a complex mix of star formation and AGN activity intermediate 
between the cold star formation dominated phase and the birth of a bright, unobscured quasar. 
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1 INTRODUCTION 

With bolometric energy outputs rivaling bright quasars, luminous 
(LIRGs) and ultraluminous (ULIRGs) infrared (IR) galaxies are 
some of the most extreme objects in the Universe (for a review, see 
ISanders & M irabel 1996). Although they make up a small fraction 
of the local IR lumin osity density, at z > 1 they become cosmo - 
logically important ihe Floc'h etal] 120051 ; iMagnelli et all 120091) . 
contributing a sig nificant fraction of the diffuse extr agalactic IR 
background (e.g.. iFixsen et alj ll99Sl ; IPei et alj|l999l) and up to 
half of the star formation in the early Universe (e.g., iBlain et all 
1999). Over the past two decades, a number of authors have ar- 
gued that these extreme systems represent an intermediate stage in 
the merger-driven tra nsformation of gas-rich quiescent disks int o 
luminous quasars (e.g. Sand ers et al. 1988a: lHopkins et all2 006c). 
As a result, a detailed understanding of (U)LIRGs - including the 
engine that powers their luminosity - is essential to any successful 
theory of galaxy formation. 

Based on a number of observational arguments, energy pro- 
duction in most (U)LIRGs appears to be dominated by star for- 
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mation (SF) both locally ( e.gjGenzel et al.lll998l;lLutz et all 1998h 
and at high redshift (e.g., I Alexander et alj|2005l ; IPope et alJboOS ; 
lYounger et alj|200l 12008 dlbl) . However, the "warm" IR colours 
seen i n some systems (S^m/Secium iJ 0.27: Ide Grijp et alj 19851 
fl987h are thought to betray a heavily obscured active galactic nu- 
cleu s (AGN) that contr ibutes significantly to the IR emission (see 
also iRisaliti et al.1l2000h . Therefore, warm ULIRGs may represent 
the fina l IR luminous stage im mediately prior to the birth of a bright 
quasar ISanders et al. 19883). 



Recent theoretical work has shown that mer gers can drive star- 
bursts consistent with those seen in (U)LIRGs (Mihos & Hernquistl 
1 1994 1 1994 iHopkins et~ai] l2006bl). Furthermo re, the radiative 
transfer calculations of Icha krabarti et al. (2007) demonstrate that 
warm IR colours are associated with increased AGN activity. How- 
ever, since the starburst and AGN activity are both produced via 
the same merger-driven inflow, it is possible that both heat the dust 
distribution. Here, we investigate the relative contribution of star 
formation and AGN activity to the production of warm (U)LIRGs. 
This work is organized as follows: in § [2] we summarize our 
methodology, in § [3] we present our results and discuss some po- 
tential interpretations, and in §|4]we conclude. 
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2 METHODOLOGY 

The hydrodynamical simul ations prese nted in this work were per- 
formed using GADGET-2 ( Sori ngell2005l) . a smoothed particle hy- 
dr odynamics (SPH) code usi ng the entropy conserving formalism 
of Springel & Hernquist] d2002l) . We include the effec ts of gas dis- 
sipation via radia tive cooling (see iKatz et al.lll99q) and a sub- 
resolution model dSpringel & Hemquistl 120031) for the multiphase 
interstellar medium (ISM) with a star formation (SF) t imescale 
tuned to match the local Kennicutt-Schmidt relation dSchmidtl 
Il959t l Kennicutt 1998). Finally, we include sink particles represent- 
ing supermassive black holes (SMBH) that undergo Eddi ngton- 
limite d Bondi-Hoyle-Littleton accretion (for a review, see Edgar 
|2004) and release thermal feedback into the interstellar medium 
i Springel et al. 2005). A number of previous studies have found 
that this model of AGN feedback, assuming a radiative efficiency 
of ~ 10% and thermal coupling efficiency of w 5%, yields 
AGN light curves that simult aneously reproduce both the observed 
quasar luminosity function (Hopkins et al. 2006^1 120081) and the 
full range of obser ved scalings between the final SMBH mass and 
galaxy properties dDi Matteo et a l. 2005; Robertson et al.l [2006b ; 
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The progen itor galaxies were co nstructed following the meth- 
ods described in Springel et al. (2005), and ar e similar to the Milky 
Way-like model in Robertson et al. (2006a): they are bulgeless 
disks embedded in a iHernquistl ( 1990) dark matter halo of mass 
M v i r = 1.4 x 10 12 Mq with concentration and spin consistent 
with cosmological N-body simulatio ns (cw = 9 and A = 0.033: 
Bullock et alj|200ll ; iBett et alj |2007) and baryonic mass fraction 
nib — 0.041. Each of the progenitor disk galaxies is given an ini- 
tial gas fraction of f g — 0.4 by mass. They are then placed on 
a zero-energy pa rabolic orbit, as motivated by cosmo logical N- 
body simulations (Benson 2005; Khochfar & Burkert 2006). In this 
work, we exam ine three different orbital configurations (see also 
ICox et al"1 l2006): orientation e (tilted prograde-prograde: 61 = 30, 
4>i = 60, 82 = —30, 4>2 = 45), h (coplanar prograde-prograde: 
0i = 0, (f>i = 0, 02 = 0, cf> 2 = 0), and g (tilted prograde- 
retrograde: 9i — 150 4>i = 0, 82 = —30, <$>\ = 45). The se are the 
same merger simulations used bv lNaravanan et all d2008l) to study 
the effects of AGN- and starburst- driven winds on molecular gas 
emission from IR luminous galaxies. 

Simulation snapshots were post-processed using a signifi- 
cantly updated version o f the polychrom atic Monte Carlo radiative 
transfer code SUNRISE dJonssonl 2006, updates will be described 
in Jonsson & Groves, in preparation). The adaptive grid resolu- 
tion was set to 5 kpc allowing for 9 levels of refinement with a 
lower limit of two gas smoothing lengths, which at the peak of 
the starburst corres ponded to ~ 30 pc. We use the R=3.1 Milky 
Way dust model of Weingartner & Draine (2001) with a constant 
Gas/DusQ ratio of 50 by mass, which is consistent with obser- 
vations of local ULIRGs dwilson et all l2008bl : Ifono et alj|2009h . 
We consider two different treatments of dust in the ISM: (1) the 
dust density in a given grid cell is determine d using the multiphase 
ISM model o f lSpringel & Hemquistl dlooH) assuming a negligible 
volume filling fraction for dense molecular gas - i.e., only dust in 
the diffuse phase contributes to the attenuation - and (2) a single- 



We haven chosen to set a constant Gas/Dust ratio, rather than the 
more physically motivated Metals/Du st, because it is dire ctly constrained 
by observations of (U)LIRGs (e.g., Iwilson et all l2008al) . but note that 
the IR colours for Gas/Dust=50 are consistent with those assuming Met- 
als/Dust=0.4 (Dwek 1998) to within 10%. 



phase ISM in which the dust is uniformly distributed in the grid 
cells - a configuration that h as been shown to maximize the atten- 
uation dWitt & Gordoni ri996). Furthermore, we treat the effects of 
dust self-absorption properly by iterating the dust emission and ab- 
sorption stages of the radiative transfer ca lculation until the dust 
temperature converges dMisselt et ai1l200ll) . Jonsson & Groves (in 
preparation) find that their implem entation recovers the b enchmark 
2-D radiative transfer solution of IPascucci et alj d2004l) to better 
than ~ few % for UV through millimeter wavelengths. Finally, we 
consider the emergent IR SED over 6 different sight lines, but note 
that at A > 20/im it is fairly isotropic (see Figure!]}. 

We assign the stellar particles SEDs from the library of 
I Groves et ID §008), which accounts for the effects of the HII re- 
gions and photodissociation regions (PDRs) - the transition layer 
between the HII region and th e surrounding molecular cloud - that 
surround young star clusters. iGroves et al. (2008) use input spec- 
tra computed using the Starburst99 (Leitherer e t al J 1 19991) stellar 
population synthesis code with a Kroupa (2001) initial mass func- 
tion (IMF). They calculate the transfer of line and continuum ra- 
diation, including the effects of both gas and dust, through the sur- 
rounding HII regi on (which evolves as a one-dimensional mass-loss 
bubble following ICastor et aljfl9751) using the MAPPINGS III code 
dGroves et alj|2004) . The HII regions absorb most of the ionizing 
radiation and are responsible for essentially all the hydrogen line 
emission and the hottest dust emission. As a subresolution treat- 
ment of dense molecular gas, a fraction / p dr - which is related 
to the cloud clea ring timescale as / p d r = exp (— t/r c i car ) (see 
IGroves et al1l2008h - of the emergent luminosity of young stellar 
clusters are surrounded by PDRs. The PDRs absorb a large fraction 
of the non-ionizing UV photons and potentially contribute signifi- 
cantly to the PAH and FIR emission; this parameter encodes the ef- 
fects of obscuration and reprocessed emission owing to cold, dense 
clouds immediately surrounding the young star-forming regions. 

Finall y, the SMBH particle s are assigned an intrinsic SED 
from the Hopk ins et ail (2007b) library of unobscured quasar 
templates according to their bolometric luminosities, Lagn = 
cMbhc 2 , where Mbh is the accretion rate a nd e = 0.1 is the ra- 
diative efficiency corresponding to a standard Shakur a~& Sunvaevl 
dl973l) thin disk. This empirical model SED includes absorption 
and emission - including the torus contribution - arising on scales 
below the resolution of the SUNRISE grid ( <i 30 pc, see above). 
It does not, however, include significant additional obscuration on 
those scales from, e.g., an intrinsically Compton-thick AGN. 



3 RESULTS AND DISCUSSION 

In FigureQjwe present the 25-to-60 /im flux density ratio typically 
used to distinguish warm (U)LIRGs, including runs both with and 
without the AGN contributio n, and in Figure [2] we sho w the total 
IR luminosity (as defined in Sanders & Mirabel 1996) as a func- 
tion of simulation time. We choose a fiducial value of / p( jr = 0.2 
- representative of dense molecular gas around young stellar clus- 
ters with a 20% covering fraction - which is roughly consistent 
with some observations of the nuclear starbursts in ULIRGs (e.g., 
iDownes & Solomon| [l998). We find that, given this modest cover- 
ing fraction, while the AGN contributes additional heating, our sim- 
ulations produce a warm (U)LIRG even when only stellar emission 
is considered. Therefore, we find that that stellar emission alone 
may be sufficient to produce warm IR colours while this system is 
a ULIRG. 

These results are consistent with broadly attributing the warm- 
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Simulotion Time [Myr] Wavelength [^m] 

Figure 1. LEFT: Time evolution of the 25-to-60 fim flux density ratio during the final stages of the merger for the e orbital orientation (see §l2land lCox et alj 
2006) assuming a 20% covering fraction of PDRs - cold, dense gas immediately surrounding the young stellar particles - and a negligible volume filling 
fraction of molecular gas in the diffuse ISM. The grey and red shaded areas show the range over six different sightlines when excluding and including the 
AGN respectively; the black solid line indicates the median value for each. The das h-dotted black line indicates results including an AGN but assuming a 
single-phase ISM, which has been shown to maximize the attenuation (Witt & Gordon 1996). Dotted lines with a star and circle identify the time of peak of 
star formation and AGN activity re spectively, and the d ashed black line shows the traditional criterion for warm IR objects 525/jm/ Sao/nn ?Z 0.27 (spectral 
index Q(25Atm, 60/im) > -1.5: lde Grijp et alj[l98lll987h . We find that given a low volume filling fraction of dense molecular gas, our simulations and 
radiative transfer calculations produce warm (U)LIRGs. RIGHT: Examples of the emergent SEDs - normalized at 25/^m- during the cold (blue) and warm 
(red) evolutionary phases including both AGN and star formation. The shaded region indicates the range over all sight-lines, and the black solid line shows 
the mean value. 



ing of the IR SED to AGN activity. However, they also demonstrate 
that warm IR sources do not necessarily contain a significant con- 
tribution from an AGN. Rather, these objects necessarily contain 
a mix of AGN and star formation that both contribute UV pho- 
tons that, after being reprocessed by the intervening ISM, lead to 
a hot dust component at 25 ;S A <J 60 fim in their SED. Figure [3] 
shows this explicitly: while there is a strong correlation between 
S^snm/SWm an£ * mtrms i c AGN luminosity (Lagn) - a rank 
correlation coefficient analysis shows that orbits e, h and g have 
highly significant (> 6a) strong positive correlations of p = 0.8, 
0.8, and 0.9 respectively - there is also significant scatter about the 
mean trend. This scatter is indicative of contributions of star for- 
mation to the warming of the IR SED. Furthermore, the 25-60 (im 
colours are produced by the large-scale ISM at all times; even if all 
of the intrinsic AGN luminosity were reprocessed on small scales 
and reemitted at 25-60 /im , it would still only contribute 10-15% 
of the emergent luminosity at those wavelengths. 

At the same time, the maxima in S l 25 /J m/S'60/jm do not oc- 
cur at the absolute peak of AGN activity. In our simulations, the 
emission at 25 and 60 /im is always dominated by dust in the dif- 
fuse ISM. Therefore, these simulated warm IR colours cannot be 
attributed directly to the intrinsic AGN emission. Unfortunately, 
the complex geometry of the obscuring material and source dis- 
tribution make it difficult to isolate the origin of this effect, and a 
number of processes may contribute. However, it may owe in part 
to a combination of dust self-absorption and gas depletion. During 
the peak of AGN activity, the gas in the nucleus is very dense and 
may be optically thick to its own emission; this will tend to lower 
the temperature despite a higher input of UV photons. At the same 
time, gas depletion due to the nuclear starburst and outflows driven 
by AGN feedback will tend to reduce the dust mass in the nucleus 
after AGN activity peaks. Since the system will still reprocess al- 
most all of the input UV/optical into the IR, and the IR luminous 
scales roue hly as Lib. ~ M d T^ +f3 (where (3 « 1.5 is the dust 
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Figure 2. Same labeling as Figure \T\ for the total infrared luminosity (8 — 
lOOO^m, as defined in Sanders & Mirabel fl996T) . also for the e orbital 
orientation. The shaded region indicates when the object satisfies the warm 
criterion of lde Grijp et1uTll985l.ri987t) . and it is considered a ULIRG at 
any time when L[r > 10 12 Lq. We find that the simulated object is 
warm during the ULIRG phase. Furthermore, we find that at times the AGN 
contributes up to ss 30 — 50% of the peak IR luminosity during this phase. 

emissivity; see e.g.. |Pe Breuck et"al]|2003l ; I Younger et alj |2008b). 
this will tend to raise the dust temperature at later times. 

We have also performed some preliminary tests in which we 
assume a single-phase ISM (see Figures[T]and[2), which h as been 
shown to maximize the attenuation ( Witt & Gordonl ll996n . In this 
scenario, the combined effects of a higher dust mass - which re- 
quires a colder effective dust temperature at fixed Lir (see above) 
- and self-absorption - when the dust is optically thick to its own 
emission, also leading to lower effective dust temperature - dramat- 
ically suppress warm IR colours both with and without an AGN 
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Figure 3. The 25-to-60 fim flux density ratio as a function of AGN 
luminosity assuming a radiative effi ciency of e = 0.1 corresponding 
to standard IShakura & Sunvaevl J 1973b thin disk accretion. Included are 
result s from three different orbital orientations (see § [5] and ICox et al] 
2006): configurations e (black circle), h (blue star), and g (red trian- 
gle). The dashed black line shows the traditional cutoff for warm IR ob- 
jects 525um/'S'60mn <D -27 (spectral index a(25fj,m, 60/jm) > —1.5: 
Ide Grijp et al1ll985lll987l) . We find a strong positive correlation between 
>525/im /SeOfim and Lagn in all cases, but with significant scatter about 
the mean trend. The significant scatter in this scaling indicates that while the 
warm phase is broadly correlated with increased AGN activity, a number of 
factors are significant including UV emission from young stellar popula- 
tions and the spatial distribution of dust. 



contribution while reaching very similar total IR luminosities; in 
this maximal attenuation case, even the AGN is insufficient to pro- 
duce warm colours. Though a more sophisticate d implementation 
of co l d dense gas in the ISM is warranted (e.g.. Narayanan et al.l 
2008; Li et al . 2008), these preliminary tests suggest that the pro- 
duction of warm IR colours is very sensitive to the effects of atten- 
uation by molecular gas along the line of sight. Furthermore, while 
assuming a single-phase ISM may seem an extreme limiting case, 
observations of local ULIRGs indicate that the covering fractions 
of optically thick molecular gas in the nucl ear starbursts powerin g 
local ULIRGs may be close to unity (e.g. JSakamoto et"ai]|l999h . 
This potentially complicates the interpretation of IR colours as an 
indicator of AGN activity. 

To test the predictions of our models observationally, we com- 
pare our simulated systems - including two additional orbital con- 
figurations - with hard X-ray (HX: 2-10 keV) ob servations of anal- 
ogous objects (see Figure l4l iRisaliti et al.ll200ob . We find that our 
simulations - assumin g a radiative efficiency of 10% and the in- 
trinsic SED library of iHopkins et alj d2007bl) for the AGN - are 
consistent with the observed scatter in Cir within a narrow range 
in Lhx /Liu. The relatively small dynamical range of Lhx /Lir 
seen in our simulations is not surprising; we have only sampled a 
single mass encounter over a small number of orbital configurations 
and parameter choices. Indeed, a thorough study of the effects of 
varying the Gas/Dust ratio, dense molecular gas along the line-of- 
sight, galaxy/halo mass, and others are all topics we plan to study in 
detail in future work. Furthermore, there could be additional obscu- 
ration near the AGN on scales below our spatial and mass resolu- 
tion or owing to dense gas along the line of sight would yield signif- 



Figure 4. Correlation between Lhx I Lir - where Lhx is the hard X- 
ray (2-10keV; estimated assu ming a radiative efficie ncy of 10% and the 
intrinsic AGN SED library of Hopkins et al. 2 007bl) and L IR is the to- 
tal (8-1000^tm) IR luminosity - and Cir = 2 X S^Sum /'S'60 um for our 
simulations compared to the observations of Risaliti et al. (2000). T he gre y 
shaded area shows the range in observed systems from Risaliti etal. (2000). 
Show n are results from three orbital configurations (see £j|2land lCox et alj 
|2006|) - e (black circle), h (blue star) and g (red triangle) - and / p dr = 0-2. 
We furthermore restrict ourselves to times when the simulated systems are 
at > 0.25 X max(SFR) or > 0.25 X max(LAGN ) - though this does not 
significantly change the region occupied by our simulations - to exclude 
times when the progenitor galaxies are widely separated; our models con- 
cern merger-driven starbursts rather than isolated disks. Fi nally, we also 
indicate the location of unobscured quasar templates from Hopkins et al. 
(2007b) for L AGN = lO 11 " 13 Z/q . We find that our merger-driven evolu- 
tionary models of warm (U)LIRGs are consistent with the observed scatter 
in 25-to-60 fim flux density ratio at fixed Lhx / Lir. 



icantly higher column densitiefl similar to those which have been 
observed locally (e.g.. iMaiolino et all 1 19981 ; iBassani et al. I [19991 
which would tend to lower both Lhx/Lir and S2511U1/ Seo^m- 
Conversely, lower column densities along particular sight lines 
could tend towards higher Lhx/Lir and S^SMm/Seo/mi- Finally, 
gas-rich isolated star -forming disks with radiatively inef ficient or 
low-level AGN (e.g.. lHopkins & Hernquis3l2006ll2008i3 lbl) could 
have low Lhx/Lir and a range of S25 M m/£Wm colours. How- 
ever, our models of a merger-driven evolutionary scenario for the 
production of warm (U)LIRGs a re consistent with th e constraints 
imposed by the observations of Ris aliti et alj d2000h and others. 
However, at the same time our simulations demonstrate that the 
correlation between Cir and Lhx/Lir should not be taken to 
show a direct link between warm colours and AGN activity; rather, 
their origin is inherently ambiguous, produced by star formation, 
AGN activity, or both in combination. 



4 CONCLUSION 

We present the results of radiative transfer calculations, including 
a self-consistent treatment of dust heating, on simulations of gas- 
rich major mergers. Our models produce warm IR SEDs in systems 
analogous to local (U)LIRGs, and we find that the the warm phase 



In our simulations, the optical dep th of the HX to photoelectric absorp- 
tion ( Morrison & McCammon 1983) is r < 0.05 at all times of interest. 
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is broadly correlated with AGN activity. Furthermore, the mid-IR 
emission arises almost entirely from the diffuse ISM rather than the 
intrinsic AGN SED. But, we show that star formation alone may be 
sufficient to produce warm IR colours, which suggests that they 
may be a more ambiguous probe of AGN activity. However, these 
results are sensitive to the multiphase structure of the ISM; when 
we c onsider a single-phas e ISM, which maximizes the obscura- 
tion iWitt&Gordor]|l996h . even a significant AGN contribution 
will not produce warm IR colours. Finally, we find that the HX and 
IR properties - and specifically the observed scatter in the 25-to- 
60 /im flux density ratio at fixed Lhx /Lir - of local (U)LIRGs 
are consistent with our models of a merger-driven evolutionary 
scenario in which both star formation and AGN activity probably 
contribute significantly to the production of warm (U)LIRGs, and 
that these objects are intermediate between the star formation dom- 
inated cold phase and the birth of a bright, unobscured quasar. 
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